﻿using System;
using System.Collections.Generic;
using Allegro.MathInterface;
using Allegro.MathLib;
using Allegro.TaskExecution;

namespace Allegro.Mathlib
{
    /// <summary>
    /// Generic equation solver for finding multiple solution dispersed over a predefined
    /// hyper cube
    /// </summary>
    public class MultiSolutionEquationSolver
    {
        private IEquationSetSolver _equationSolver;
        private double[] _lowerBounds;
        private double[] _upperBounds;
        private LatinHypercubeSampler _hypercubeSampler;
        private VectorMap _function;
        private double[] _steps;
        private List<double[]> _solutions;
        private int _maxTrials;
        private int _maxSolutions = int.MaxValue;
        private int _solutionCount = 0;

        /// <summary>
        /// 
        /// </summary>
        /// <param name="f"></param>
        /// <param name="equationSolver">Solver used</param>
        /// <param name="trialSteps"></param>
        /// <param name="lowerBounds">Lower bounds, one per variable</param>
        /// <param name="upperBounds">Upper bounds, one per variable</param>
        /// <param name="nSamples">Number of samples per variable</param>
        public MultiSolutionEquationSolver(VectorMap f, IEquationSetSolver equationSolver, double[] trialSteps,
            double[] lowerBounds, double[] upperBounds, int nSamples)
        {
            if (nSamples < 2)
            {
                string msg = string.Format("MultiSolutionEquationSolver: Number of samples per variable must be >= 2");
                Log.Error(msg, null, LogEventTypes.OtherEvent);
                throw new ArgumentException(msg);
            }
            _function = f;
            _steps = trialSteps;
            _equationSolver = equationSolver;
            _lowerBounds = lowerBounds;
            _upperBounds = upperBounds;

            int nParameters = lowerBounds.Length;
            _maxTrials = nSamples * nParameters;
            double[][] samples = new double[nParameters][];
            for (int i = 0; i < nParameters; i++)
            {
                samples[i] = new double[nSamples];
                for (int j = 0; j < nSamples; j++)
                {
                    samples[i][j] = lowerBounds[i] +  (upperBounds[i] - lowerBounds[i]) * ((double)j) / ((double)(nSamples - 1));
                }
            }
            _hypercubeSampler = new LatinHypercubeSampler(nParameters, nSamples, samples);

            _solutions = new List<double[]>();
        }
        /// <summary>
        /// Maximum number of solutions sought. Default is indefinite
        /// </summary>
        public int MaxSolutions
        {
            get { return _maxSolutions; }
            set{ _maxSolutions = value;}
        }
        /// <summary>
        /// Find next solution. If false is returned, no solution was found 
        /// </summary>
        /// <param name="root"></param>
        public bool NextSolution(ref double[] root)
        {
            double[] x = null;
            bool rootFound = false;
            int iTrial = 0; // Trial counter
            while (_solutionCount < _maxSolutions 
                && iTrial < _maxTrials 
                && _hypercubeSampler.More 
                && !rootFound)
            {
                x = _hypercubeSampler.NextSample();
                iTrial++;
                try
                {
                    _equationSolver.FindRoot(_function, x, _steps, ref root);
                    if (_RootAlreadyFound(x))
                    {
                        rootFound = false;
                        break;
                    }
                    rootFound = true;
                    _solutionCount++;
                }
                catch (Exception ex)
                {
                    rootFound = false;
                }              
            }
            
            return rootFound;
        }

        private bool _RootAlreadyFound(double[] xValue)
        {
            bool found = false;
            int n = xValue.Length;
            foreach(double[] root in _solutions)
            {
                found = true;
                for (int i = 0; i < n; i++)
                {
                    if (Math.Abs(xValue[i] -root[i]) > _equationSolver.ToleranceX)
                    {
                        found = false;
                        break;
                    }
                }
                if(found)
                    break;
            }
            if(!found)
                _solutions.Add(xValue);
            return found;
        }
    }
}
